Cite this Article
Wang Zhan-Feng, Zhang Zhi-Qiang, Fu Yi-Ben, Wang Peng-Ye, Xie Ping. Brownian ratchet mechanism of translocation in T7 RNA polymerase facilitated by a post-translocation energy bias arising from the conformational change of the enzyme. Chinese Physics B, 2017, 26(3): 030201  
Permissions
Brownian ratchet mechanism of translocation in T7 RNA polymerase facilitated by a post-translocation energy bias arising from the conformational change of the enzyme
Wang Zhan-Feng, Zhang Zhi-Qiang, Fu Yi-Ben, Wang Peng-Ye, Xie Ping †
Key Laboratory of Soft Matter Physics and Beijing National Laboratory for Condensed Matter Physics, Institute of Physics, Chinese Academy of Sciences, Beijing 100190, China

 

† Corresponding author. E-mail: pxie@aphy.iphy.ac.cn

Abstract

T7 RNA polymerase can transcribe DNA to RNA by translocating along the DNA. Structural studies suggest that the pivoting rotation of the O helix in the fingers domain may drive the movement of the O helix C-terminal Tyr639 from pre- to post-translocation positions. In a series of all-atom molecular dynamics simulations, we show that the movement of Tyr639 is not tightly coupled to the rotation of the O helix, and that the two processes are only weakly dependent on each other. We also show that the internal potential of the enzyme itself generates a small difference in free energy (ΔE) between the post- and pre-translocation positions of Tyr639. The calculated value of ΔE is consistent with that obtained from single-molecule experimental data. These findings lend support to a model in which the translocation takes place via a Brownian ratchet mechanism, with the small free energy bias ΔE arising from the conformational change of the enzyme itself.

PACS: 02.70.Ns;05.40.-a;05.40.Jc;87.15.kj
Keyword:RNA polymerase;molecular dynamics simulation;molecular motor;Brownian ratchet
1. Introduction

RNA polymerases (RNAPs) are enzymes that transcribe genes from DNA onto strands of RNA. Two families of these enzymes exist in nature: single-subunit RNAPs, found in mitochondria, some chloroplasts, and phages, and multisubunit versions, found in archaea, bacteria, animal viruses, eukaryotic, and plant nuclei.[1] The transcription elongation by RNAP is a processive process: one enzyme can catalyze thousands of nucleotide addition reactions without dissociating from the DNA/RNA substrate.[2–4] During this processive transcription elongation, the translocation of RNAP along the DNA template is “powered” by the chemical energy of nucleotide triphosphate (NTP) hydrolysis. The manner in which RNAP converts this chemical energy into mechanical translocation and the molecular mechanisms underlying this translocation have been the subject of intense investigation.

The molecular architecture of the single-subunit T7 RNAP resembles the right-hand-like configuration of DNA polymerase I, with the active site for nucleotide incorporation located on the palm of the hand.[5–7] An α-helix (called O helix) in the fingers domain, which abuts the active site, makes pivoting rotation as the fingers domain itself transitions between open and closed conformations during each nucleotide addition or elongation cycle.[8] A series of structural studies on T7 RNAP supported the existence of a power-stroke mechanism of translocation.[8–10] In this mechanism, the release of pyrophosphate (PPi) triggers the pivoting rotation of the O helix along with the outward rotation of the fingers domain, causing the residue Tyr639 on the Cterminal of the O helix to move towards the primer terminus of the hetero-duplex product and thus facilitating the translocation of the DNA-RNA via the interaction between the side chain of Tyr639 and the nascent hybrid base pair. In contrast, both biochemical[11–14] and single-molecule[15–17] experiments have indicated that thermal fluctuations play active roles in the translocation, where the translocation is driven thermally and the substrate binding is considered to act as the “pawl”, biasing the rapid forward-backward translocation equilibrium between the pre- and post-translocation states. Interestingly, the single-molecule data also revealed a small free energy bias toward the post-translocation state over the pre-translocation state.[17,18] To resolve the conflicting views between the structural data[8] and the single-molecule experimental data,[17,18] a hybrid model was recently proposed, in which the facilitated translocation or power stroke coexists with the Brownian-ratchet-driven motions.[19] However, further structural studies are still needed to test the model. Moreover, the origin of the small difference in free energy between the pre- and post-translocation states remains to be clarified.

The static structures of the elongation complex in pre- and post-translocation states determined using x-ray crystallography have formed the basis of our understanding of the translocation mechanism.[8–10] For a more complete understanding of this process, it is crucial to determine the dynamic structures that form during the transition from pre- to post-translocation states, thus revealing how different components involved in the translocation coordinate with each other to induce or facilitate the translocation. This work explores the dynamics of T7 RNAP by means of all-atom molecular dynamics (MD) simulations, which are used extensively to study the dynamics of RNAPs[20–24] and DNA polymerases.[25–27] A novel method is applied to prevent the translocation of atomistic MD simulations from taking too long (on the millisecond timescale or even longer). This simulates the movement of Tyr639, the pivoting rotation of the O helix, and the outward rotation of the fingers domain during the transition from the pre- to post-translocation states, thus showing how these components are coordinated with each other and what role each component plays in the translocation process.

Unexpectedly, the results show that the movement of Tyr639 is only weakly dependent on the outward rotation of the O helix and the fingers domain, suggesting that the motion of Tyr639 is not driven by the pivoting rotation of the O helix. The results also show that the internal potential of the enzyme itself, which is dependent only on the relative distance between Tyr639 and the palm domain, produces a small free energy bias toward the post-translocation state rather than toward the pre-translocation state. Since the interactions between Tyr639 and primer terminus of the hybrid in pre- and post-translocation states are nearly the same (see below), the translocation could be facilitated mainly by the internal potential of the enzyme itself. Thus, the translocation could be facilitated mainly by the internal potential of the enzyme itself. The results also show that the calculated difference between the internal potential energy at pre-translocation state and that at post-translocation states is quantitatively supported by the available single-molecule experimental data. Moreover, an intermediate state of translocation is revealed in which Tyr639 is located in between the pre- and post-translocation positions. These findings may significantly advance the understanding of the translocation mechanism.

2. Results

The representative structures in an elongation cycle, as reported in available studies, are shown in Fig. 1.[8–10] In the pre-translocation complex with Tyr639 away from the active site and the fingers domain in the closed conformation (Fig. 1(b)), the strong interaction between the side chain of Tyr639 and the (i + 1)-th base pair at the primer terminus keeps Tyr639 in a fixed position relative to the DNA-RNA hybrid. Consequently, only the palm and thumb domains can move relative to the DNA-RNA toward Tyr639 and the fingers domain, becoming the post-translocation state (Fig. 1(c)). Due to the strong interaction between the palm and thumb domains and the DNA-RNA, the movement of the palm and thumb domains relative to the DNA-RNA is estimated to require a millisecond or even longer,[28] a timescale too long for atomistic MD simulations of this motion. However, if the (i + 1)-th base pair and PPi (Fig. 1(b)) are removed, the movement of Tyr639 towards the palm domain accelerates significantly, and can be simulated easily using the all-atom MD method. By doing so, the movement of Tyr639 and the rotation of the O helix along with the rotation of the fingers domain can be simulated, which allows the study of the dynamic relationship between them during the translocation process. After Tyr639 moves into position with its side chain at the active site and the conformation of the fingers domain changes from closed to open, the elongation complex enters the post-translocation state (Fig. 1(c)) except that the side chain of Tyr639 forms interaction with the i-th base pair.

Fig. 1. (color online) T7 RNAP translocation-elongation cycle. (a) A substrate insertion complex (PDB 1S76) with the fingers domain in closed conformation and residue Tyr639 away from the active site (referred to as PRE position). The template nucleotide to be transcribed is numbered i + 1. The NTP binding site (active site) is shown in a red oval. (b) A pre-translocation complex (PDB 1S77) with the fingers domain in closed conformation and residue Tyr639 in PRE position. (c) A post-translocation complex (PDB 1H38) with the fingers domain in open conformation and residue Tyr639 close to the active site (referred to as POST position). (d) A pre-insertion complex (PDB 1S0V) with the fingers domain in open conformation and residue Tyr639 in POST position.

As shown in the simulations, when there is no interaction with the DNA-RNA base pair, the side chain of Tyr639 is highly mobile, but the CA atom of Tyr639 still has definite positions as in the case in the presence of the interaction. For this reason, the position of Tyr639 is determined using that of its CA atom. Tyr639 is therefore defined to be in the “PRE” position when its CA atom is in the pre-translocation conformation, away from the active site along the translocation direction (Fig. 1(b)), and in the “POST” position when its CA atom is in the post-translocation conformation, close to the active site (Fig. 1(c)). The movement of Tyr639 was characterized using the change in position of the CA atom of Tyr639 along the line connecting it in PRE position (PDB 1S76 or 1S77) and POST position (PDB 1H38), after superimposing the respective palm domains with UCSF Chimera package.[29,30] This change in the position of the CA atom of Tyr639 is denoted by a distance . From the superimposed structures, it was determined that the transition from the pre-translocation conformation with PRE Tyr639 to post-translocation conformation with POST Tyr639 corresponds to nm. Thus, in our simulations, when Tyr639 moved to positions with nm and then stayed there for a time period τ ≥ 6 ps, it was considered that Tyr639 has moved to the POST position. Note that taking other values of period τ has no effect on our studies of the dynamic relationship between the movement of Tyr639 and the rotation of the O helix. As indicated by the structural data,[8] the restriction of the O helix in the closed conformation before the release of PPi was caused primarily by the interaction of PPi with Arg627 located at the N-terminus of the O helix. For this reason, the rotation of the O helix after the release of PPi was characterized by the change in the position of CA atom of Arg627 along the line that connects the CA atom of Arg627 in pre-translocation conformation (PDB 1S76 or 1S77) and that in post-translocation conformation (PDB 1H38), after superposition of their respective palm domains. This change in the position of the CA atom of Arg627 is here denoted by a distance . Note that using distance to characterize the overall rotation of the O helix is equivalent to using the rotation angle of the O helix. From the superimposed structures, it was found that the rotation of the O helix and fingers domain from the closed to open conformation corresponds to nm. We define the O helix and fingers domain as being in the open conformation when the CA atom of Arg627 has moved by nm and then stayed there for a period τ ≥ 6 ps. As shown by the structural data,[8] the characteristic difference between the pre-translocation and post-translocation conformations of the enzyme is represented by the difference in the distance between Tyr639 together with the fingers domain and other domains along the translocation direction (comparing Fig. 1(b) with Fig. 1(c)). Thus, throughout our simulations some residues of palm domain were restrained to prevent the overall drift and rotation of the complex (see simulation methods).

2.1. Intermediate conformations revealed for Tyr639 and the fingers domain

The simulation was started with the pre-translocation complex (Fig. 1(b), PDB 1S77) by removing PPi and converting 3′-deoxyadenosine to adenosine, where initially the fingers domain is in the closed conformation and Tyr639 is in the PRE position (called system 1). As described above, due to the strong interaction between the side chain of Tyr639 and the DNA-RNA base pair at the primer terminus, Tyr639 remained fixed relative to the DNA-RNA, and due to the strong interaction between the palm and thumb domains with the DNA-RNA, the movement of the DNA-RNA relative to the palm and thumb domains was estimated to be on the millisecond timescale or even longer. For this reason, in the current simulations within a 10 ns period, it is expected that neither the DNA-RNA nor Tyr639 would move relative to the palm and thumb domains. These simulations showed that, within the 10 ns period, Tyr639 remained stable in the PRE position and the position of 3′-end of RNA also remained unchanged (Fig. 2(a)). Interestingly, Arg627 or O helix moved rapidly to an intermediate position (Fig. 2(a)). This position is here called the “semi-closed” conformation of the O helix and fingers domain. The structure is shown in Fig. 2(b). A statistical analysis of the position distribution of Tyr639 (Fig. 2(c)) and Arg627 (Fig. 2(d)) was performed with 20 trajectories with a simulation period of 10 ns. As shown in Fig. 2(c), Tyr639 remained at the PRE position and its fluctuation around the PRE position was only about ±0.1 nm. This is inconsistent with the hybrid model, which argues that the 3′-end of the RNA is quite flexible and only a small amount of energy is required to move Tyr639 from PRE to POST positions in the pre-translocation complex.[19] As shown in Fig. 2(d), Arg627 in the semi-closed position deviated from the closed position by about 0.05 nm. The existence of this semi-closed conformation may be explained by the observation that after the release of the PPi, the O helix and the fingers domain tend to rotate slightly towards the open conformation.

Fig. 2. (color online) Intermediate position of Tyr639 and an intermediate semi-closed conformation of the fingers domain. (a) The time series of the change in position of Tyr639 relative to the PRE position (red line) and that of Arg627 relative to the position when the O helix and fingers domain are in the closed conformation (blue line) for system 1 where PPi is removed and 3′-deoxyadenosine is converted to adenosine in the pre-translocation complex (PDB 1S77). The latter trace is offset by −0.2 nm on the y axis for clarity. (b) A snapshot of the structure produced from the trajectory shown in panel (a), where Tyr639 is in PRE position and the O helix and fingers domain are in the semi-closed conformation (yellow). For comparison, the structure of the enzyme in the pre-translocation complex with Tyr639 in the PRE position and the O helix and fingers domain in the closed conformation (red) and the structure of the enzyme in the post-translocation complex with Tyr639 in POST position and the O helix and fingers domain in the open conformation (green) are also shown. (c) The position distribution of Tyr639 for system 1 with 20 trajectories with a simulation period of 10 ns. The blue line represents the Gaussian fit. (d) The position distribution of Arg627 for system 1 with 20 trajectories for which the simulation period is 10 ns. The blue line represents the Gaussian fit. (e) The time series of the change in position of Tyr639 relative to the POST position (red line) and that of Arg627 relative to the position when the O helix and fingers domain are in the open conformation (blue line) for system 2 which is the post-translocation complex (PDB 1H38). The latter trace is offset by −0.2 nm on the y axis for clarity. (f) A snapshot of the structure obtained from the trajectory shown in panel (e), where Tyr639 is in the ITM position and the O helix and fingers domain are in the open conformation (yellow). The structures of the enzyme shown in red and green are the same as those in panel (b). (g) The position distribution of Tyr639 for system 2 with 20 trajectories for which the simulation period is 10 ns. Yellow and green lines represent Gaussian fits; the blue line represents the sum of two Gaussians. (h) The position distribution of Arg627 in system 2 with 20 trajectories with a simulation period of 10 ns. The blue line represents the Gaussian fit.

The post-translocation complex was also simulated (Fig. 1(c), PDB 1H38). Initially, the fingers domain was in the open conformation and Tyr639 was in the POST position (called system 2). The simulations showed that the interaction between the side chain of Tyr639 and the base pair at the primer terminus was strong enough to prevent Tyr639 from moving further toward the primer terminus (Fig. 2(e)). As shown in the time series of the motion of Tyr639 (Fig. 2(e)), in addition to the POST state, there was also an intermediate (ITM) state where Tyr639 was located between the POST and PRE positions. The structure of the ITM state is shown in Fig. 2(f). The statistical results regarding the position distribution of Tyr639 with 20 trajectories or realizations within the 10 ns simulation period indicated the existence of the ITM state (Fig. 2(g)). Statistical analysis was also performed on the position distribution of Arg627 (Fig. 2(h)), showing that the open conformations of the O helix and fingers domain are stable.

Taken together, these simulations showed that there exists an ITM state in which Tyr639 is located between the PRE and POST positions and that there also exists a semi-closed conformation of the O helix and fingers domain after the release of PPi, with the O helix and fingers domain being slightly more open than in the closed conformation.

2.2. Weak coupling between the PRE to POST movement of Tyr639 and the outward rotation of the O helix and fingers domain

The time series of the motion of Tyr639 and that of O helix show that the motion of Tyr639 could not be tightly coupled to the rotation of the O helix and the fingers domain (Figs. 2(a) and 2(e)). To test this, simulations incorporating an insertion complex were performed (Fig. 1(a), PDB 1S76) by removing the α-methylene ATP and the template nucleotide i + 1 to allow the movement of Tyr639 from PRE to POST, where the fingers domain is initially in the closed conformation and Tyr639 is initially in the PRE position (called system 3). The simulations clearly showed that the movement of Tyr639 to the POST position does not occur simultaneously with the outward rotation of the O helix and fingers domain to the open conformation, as shown in Fig. 3(a). While Tyr639 moved to the POST position, the O helix had not yet rotated to the open conformation. The final structure with Tyr639 in the POST position and the O helix and fingers domain in the semi-closed conformation is shown in Fig. 3(b). Out of 20 realizations, all of which had a simulation period of 10 ns, Tyr639 moved to the POST position and then stayed there for a long time in 19. In only 1 realization did Tyr639 not move to the POST position within 10 ns (Table 1, system 3). During the motion of Tyr639 from the PRE to the POST position, there was an ITM state for which Tyr639 was located about 0.2 nm away from the PRE position (Fig. 3(a)). The existence of the ITM state was also clearly visible in the statistical analysis of the position distribution of Tyr639 for all 20 realizations (Fig. 3(c)). Data from the simulations was used to calculate the mean time for Tyr639 to move from the PRE to POST positions. This was found to be about 1628 ps (supplemental information, text S1).

Fig. 3. (color online) The movement of Tyr639 and the rotation of the O helix are weakly coupled to each other. (a) The time series of the change in position of Tyr639 relative to the PRE position (red line) and that of Arg627 relative to the position when the O helix and fingers domain are in the closed conformation (blue line) for system 3 where the αβ methylene ATP and the template nucleotide i + 1 are removed in the insertion complex (PDB 1S76). The latter trace is offset by −0.2 nm on the y axis for clarity. (b) A snapshot of the structure obtained from the trajectory shown in panel (a), where Tyr639 is in POST position and the O helix and fingers domain are in the semi-closed conformation (yellow). The structures of the enzyme shown in red and green are the same as those in Fig. 2(b). (c) The position distribution of Tyr639 for system 3 with 20 trajectories with a simulation period of 10 ns. Yellow and green lines represent the Gaussian fits; the blue line represents the sum of two Gaussians. (d) The position distribution of Arg627 for system 3 with 20 trajectories with a simulation period of 10 ns. The blue line represents the Gaussian fit. (e) The time series of the change in position of Tyr639 relative to the PRE position (red line) for system 4 where residue Arg627 was restrained in system 3. (f) A snapshot of the structure obtained from the trajectory shown in panel (e), where Tyr639 is in the POST position and the O helix and fingers domain are in the closed conformation (yellow). The structures of the enzyme shown in red and green are the same as those in Fig. 2(b). (g) The position distribution of Tyr639 for system 4 with 20 trajectories with a simulation period of 10 ns. Yellow and green lines represent Gaussian fits; the blue line represents the sum of two Gaussians.

By contrast, for all 20 realizations, neither the O helix (Arg627) nor the fingers domain had rotated to the open conformation within 10 ns, indicating that the mean time required for the outward rotation of the O helix (Arg627) and the fingers domain to the open conformation was longer than 10 ns, possibly much longer. Rather, the O helix (Arg627) and the fingers domain stayed stably in the intermediate semi-closed conformation (Fig. 3(a)). Figure 3(d) shows the statistical results of the position distribution of Arg627 with all 20 trajectories and a simulation period of 10 ns. Only one peak (corresponding to the semi-closed conformation) existed, indicating that the O helix had rotated to the semi-closed conformation very early and then stayed there for a time longer or much longer than 10 ns. These results indicated that the movement of Tyr639 from PRE to POST position is not tightly coupled to the outward rotation of the O helix and fingers domain.

Because the movement of Tyr639 to the POST position took place before the outward rotation of the O helix and fingers domain to the open conformation, the PRE to POST movement of Tyr639 could not have been caused by the outward rotation of the O helix and the fingers domain, which suggests that the power-stroke mechanism that the opening of the fingers domain together with the pivoting rotation of the O helix triggered by PPi release drive the translocation may be wrong. To further test this, in system 3, residue Arg627, which interacts with PPi, was restrained (called system 4). This mimics the case in which PPi is not released. The results of the simulation showed that the PRE to POST movement of Tyr639 can still occur without the pivoting rotation of the O helix triggered by the release of PPi (Fig. 3(e)). The structure with Tyr639 moved to the POST position was shown in Fig. 3(f). Out of 20 realizations, Tyr639 moved to the POST position and stayed there for a long time in 14 (Table 1, system 4). The statistical results of the position distribution of Tyr639 also showed the existence of an ITM state (Fig. 3(g)).

Table 1.

Time required for Tyr639 (Arg627) to move from the PRE (closed) to the POST (open) position in different simulation scenarios.

.

Taken together, the current simulations convincingly indicated that the movement of Tyr639 from the PRE to POST position is not tightly coupled to the outward rotation of the O helix and fingers domain. In addition, the existence of an ITM conformation was revealed, in which Tyr639 is located about halfway between the PRE and POST conformations. During the outward rotation of the O helix and the fingers domain from closed to open conformation, there exist corresponding intermediate semi-closed conformations for which the O helix and the fingers domain are rotated slightly outward relative to the closed conformation.

2.3. The PRE to POST movement of Tyr639 and the outward rotation of the O helix and fingers domain are weakly dependent upon each other

As seen above, without the pivoting rotation of the O helix, Tyr639 can still move from the PRE to POST positions (Figs. 3(e)– 3(g) and Table 1 for system 4). The mean time required for Tyr639 to move from the PRE to POST positions was calculated and found to be about 3017 ps when Arg627 was restrained (for system 4, supplemental information, text S1), a mean time of about 1628 ps was required for Tyr639 to move to the POST position when Arg627 was not restrained (for system 3). The effect of restraining Arg627 was found to be equivalent to a decrease in free energy of facilitating the movement of Tyr639 from PRE to POST positions by about . By fitting the available single-molecule experimental data to the Brownian ratchet model of translocation, a free energy bias of about toward the post-translocation state relative to the pre-translocation state was found (see supplemental information, text S2).[17,18] The pivoting rotation of the O helix, which is triggered by the release of the PPi, only makes a minor (about 26.8%) contribution to facilitating the forward translocation. The weak effect of restraining Arg627 on the movement of Tyr639 can be explained as follows. As shown in Fig. 3(a), after Arg627 became unrestrained, the O helix quickly rotated to the semi-closed conformation, in which the O helix is only slightly rotated relative to that in the closed conformation. This took place before Tyr639 moved from the PRE to the POST position. The O helix remained in the semi-closed conformation for much longer than the time required for Tyr639 to move from the PRE to the POST position. This means that only the rotation of the O helix from the closed to semi-closed conformation could affect the motion of Tyr639 and that the rotation of the O helix from the semi-closed to open conformation had no effect. As in the semi-closed conformation the O helix was only slightly rotated relative to the closed conformation, it would be expected that preventing the O helix from changing to the semi-closed conformation by restraining Arg627 has only a weak effect on the movement of Tyr639.

To further test the effect of the pivoting rotation of the O helix and the outward rotation of the fingers domain on the motion of Tyr639, we applied a constant force of 100 pN to seven residues (residues 654–660) in the C-terminus of the Y helix to facilitate the outward rotation of the fingers domain (called system 5). Except for the applied force, system 5 is identical to system 3. Figure 4(a) shows an example of the time series of the rotation of the O helix and the motion of Tyr639. As expected, under the applied force, the outward rotation of the O helix to the open conformation proceeded much faster (Table 1 for system 5). The simulated structure with Tyr639 moved to the POST position and Arg627 rotated to the open conformation is shown in Fig. 4(b). Of 20 realizations with the simulation time of 10 ns, the O helix rotated to the open conformation and then stayed there for a long time in 18 (Table 1, system 5). The mean time required for the O helix to rotate to the open conformation was about 4548 ps (supplemental information, text S1). Because the mean time of rotation to the open conformation with no applied force exceeds 10 ns, it was here estimated that the effect of the applied force on O helix rotation is equivalent to increasing the free energy of facilitating the rotation of the O helix from the closed to the open conformation by more than . By comparison, for all 20 realizations, Tyr639 moved to the POST position and then stayed there for a long time (Table 1, system 5). The mean time required for Tyr639 to move from the PRE to POST position was calculated to be about 1262 ps (supplemental information, text S1) under the applied force and a mean of about 1628 ps in the absence of applied force. It was here estimated that the effect of the applied force on the movement of Tyr639 movement is equivalent to increasing the free energy of facilitating movement of Tyr639 from the PRE to the POST position by about , which is much smaller than the effect on O helix rotation (with the free energy of facilitating outward rotation being increased by a value larger or much larger than 0.788kBT). This supports the hypothesis that the outward rotation of the O helix and the fingers domain has only a minor effect on facilitating the PRE to POST movement of Tyr639.

Fig. 4. (color online) The PRE to POST movement of Tyr639 and the outward rotation of the O helix and fingers domain are weakly dependent upon each other. (a) The time series of the change in position of Tyr639 relative to the PRE position (red line) and that of Arg627 relative to the position when the O helix and fingers domain are in the closed conformation (blue line) in system 5, in which a constant force of 100 pN is applied to the C-terminus of the Y helix to facilitate the outward rotation of the fingers domain in system 3. The latter trace is offset by −0.2 nm on the y axis for clarity. (b) A snapshot of the structure obtained from the trajectory shown in panel (a), where Tyr639 is in the POST position and the O helix and fingers domain are in the open conformation (yellow). The structures of the enzyme shown in red and green are the same as those in Fig. 2(b). (c) The change in position of Arg627 over time relative to the position when the O helix and fingers domain are in the closed conformation (blue line) for system 6 where Tyr639 is restrained in system 3. (d) A snapshot of the structure obtained from the trajectory shown in panel (c), where Tyr639 is in the POST position and the O helix and fingers domain are in the semi-closed conformation (yellow). The structures of the enzyme shown in red and green are the same as those in Fig. 2(b). (e) The position distribution of Arg627 for system 6 with 20 trajectories with a simulation period of 10 ns. The blue line represents the Gaussian fit.

To see if the PRE to POST movement of Tyr639 affects the outward rotation of the O helix and fingers domain, Tyr639 was restrained in system 3 (called system 6). When the residue Tyr639 was restrained, Arg627 and the O helix also rotated rapidly to the semi-closed conformation and stayed there (Fig. 4(c)), as shown in Fig. 4(d). The statistical results of the position distribution of Arg627 with 20 trajectories and a simulation period of 10 ns are shown in Fig. 4(e). The position distribution of Arg627 when Tyr639 is restrained was compared to that of Arg627 when Tyr639 is not restrained (Fig. 3(d), system 3; Fig. 4(e) for system 6). The two distributions were very similar, with maximum points at nearly the same position and similar half widths, indicating that the outward rotation of the fingers domain is only weakly dependent on the PRE to POST movement of Tyr639.

Taken together, these results indicated that the outward rotation of the fingers domain had only a minor effect on the PRE to POST movement of Tyr639 and the PRE to POST movement of Tyr639 had only a minor effect on the outward rotation of the fingers domain.

2.4. The internal potential of T7 RNAP contributes to the translocation process

As shown above, the pivoting rotation of the O helix has only a minor effect on the PRE to POST movement of Tyr639. It is not clear what force facilitates translocation in systems from which the αβ methylene ATP and the template nucleotide i + 1 have been removed. First, it is clear that the interaction between the side chain of Tyr639 and the i-th base pair can facilitate translocation. However, the energy results extracted from the energy file produced with the MD simulation show that the interactions between Tyr639 and primer terminus of the hybrid in pre and post-translocation states are nearly the same, with −7.5 kcal/mol and −7.7 kcal/mol, respectively. Thus, the interaction between the side chain of Tyr639 and the (i + 1)-th base pair remains nearly unchanged during the translocation in real situations, and consequently there must be some other effect that facilitates the PRE to POST movement of Tyr639. To determine whether the internal potential of the enzyme itself contributes to the facilitation of the PRE to POST movement of Tyr639, the DNA-RNA and αβ-methylene ATP were removed from the insertion complex (PDB 1S76, Fig. 1(a)), retaining only the enzyme (called system 7). In this system, the fingers domain was initially in the closed conformation and Tyr639 was in the PRE position. The simulations indicated that, in the absence of the interaction with the DNA-RNA base pair, Tyr639 can also move to the POST position (Figs. 5(a) and 5(b), and Table 1 for system 7). Of 20 realizations, all of which had a simulation period of 20 ns, Tyr639 moved to the POST position and then stayed there for a long time in 16. The mean time for Tyr639 to move from the PRE to POST positions was about 3948 ps (see supplemental information, text S1), which was longer than the mean of about 1628 ps in the presence of the interaction between the side chain of Tyr639 and the i-th base pair (see above). This indicated that, in addition to the interaction, there is also some other effect that facilitates the movement of Tyr639 to the POST position when the DNA-RNA is not present. Because system 7 contained only the enzyme and because the pivoting rotation of the O helix had only a minor effect on Tyr639 movement, it can be assumed that the other effects were caused by the internal potential of the enzyme itself, which is dependent solely on the relative distance between Tyr639 and the palm domain, for example, the elastic potential of the residues connecting the fingers domain containing Tyr639 to the palm domain. As shown in Fig. 5(a), in the absence of the interaction between the side chain of Tyr639 and the i-th base pair, there is also the ITM state in which Tyr639 is located about 0.2 nm from the PRE position. The existence of the ITM state can also be clearly seen in the statistical results describing the position distribution of Tyr639 (Fig. 5(c)) in which all 20 trajectories have a simulation period of 20 ns. Moreover, the statistical results on the position distribution of Arg627 (Fig. 5(d)), in which all 20 trajectories have a simulation period of 20 ns, indicated that the O helix rotated rapidly to the semi-closed conformation and then stayed there for a long time. These results are similar to those shown in Figs. 3(d) and 4(e).

Fig. 5. (color online) The translocation is facilitated by an internal potential of the enzyme itself. (a) The time series of the change in position of Tyr639 relative to the PRE position (red line) and that of Arg627 relative to the position when the O helix and fingers domain are in the closed conformation (blue line) for system 7 where the DNA-RNA and αβ methylene ATP are removed, i.e., keeping only the enzyme in the insertion complex (PDB 1S76). The latter trace is offset by −0.2 nm on the y axis for clarity. (b) A snapshot of the structure produced from the trajectory shown in panel (a), where Tyr639 is in the POST position and the O helix and fingers domain are in the semi-closed conformation (yellow). Note that because there is no interaction with the DNA-RNA base pair, the side chain of Tyr639 fluctuates considerably. However, the CA atom of Tyr639 remains stable in the POST position. The structures of the enzyme shown in red and green are the same as those in Fig. 2(b). (c) The position distribution of Tyr639 for system 7, which involves 20 trajectories with simulation periods of 20 ns each. Yellow, green, and blue lines represent Gaussian fits; the cyan line represents the sum of three Gaussians. (d) The position distribution of Arg627 for system 7 with 20 trajectories with a simulation period of 20 ns each. The blue line represents the Gaussian fit. (e) The time series of the change in position of Tyr639 relative to the PRE position (red line) for system 8 where residue Arg627 is restrained in system 7. (f) A snapshot of the structure obtained from the trajectory shown in panel (e), where Tyr639 is in the POST position and the O helix and fingers domain are in the closed conformation (yellow). Note that because there is no interaction with the DNA-RNA base pair, the side chain of Tyr639 fluctuates considerably. However, the CA atom of Tyr639 remains stable in the POST position. The structures of the enzyme shown in red and green are the same as those in Fig. 2(b). (g) The position distribution of Tyr639 for system 8 with 20 trajectories with a simulation period of 20 ns each. Yellow, green, and blue lines represent Gaussian fits; the cyan line represents the sum of three Gaussians.

Similar to the studies for systems 3 and 4 where the DNA-RNA hybrid was present, in system 7 residue Arg627 was also restrained (called system 8) to facilitate direct observation of the effect of the pivoting rotation of the O helix, which is triggered by PPi release, on Tyr639 movement. As expected, Tyr639 can also move to the POST position (Figs. 5(e) and 5(f)). Of 20 realizations with a simulation period of 20 ns, Tyr639 moved to the POST position and stayed there for a long time in 16 (Table 1, system 8). The mean time required for Tyr639 to move to the POST position was calculated to be about 7452 ps (supplemental information, text S1). The results of systems 7 and 8 revealed that the effects of restraining Arg627 on the movement of Tyr639 are equivalent to decreasing the free energy of facilitating movement of Tyr639 from PRE to POST positions by about , which is similar to the estimated using the results of systems 3 and 4, in which the DNA-RNA hybrid was present (see above). The statistics for the position distribution of Tyr639 for system 8 (Fig. 5(g)) were similar to those for system 7 (Fig. 5(c)). These results further supported the conclusion that the pivoting rotation of the O helix had only a minor effect on facilitating the PRE to POST movement of Tyr639. Rather it was facilitated mainly by the internal potential of the enzyme itself, which was dependent only on the relative distance between Tyr639 and the palm domain.

Taken together, the results shown in this section indicated that, for the systems from which the αβ methylene ATP and the template nucleotide i + 1 had been removed, regardless of the interaction between the side chain of Tyr639 and the i-th base pair, the internal potential of the enzyme itself—the elastic potential of the residues connecting the fingers domain containing Tyr639 to the palm domain — also contributes to the movement of Tyr639 from the PRE to the POST position. It is here noted that in real situations in which the nascent (i + 1)-th base pair is present (Fig. 1), during the transition from the pre- to post-translocation states (Figs. 1(b) and 1(c)), the interaction between the side chain of Tyr639 and the nascent (i + 1)-th base pair remains nearly unchanged. In this way, the translocation in real situations is facilitated primarily by the internal potential of the enzyme.

2.5. Agreement of calculated T7 RNAP internal potential with data derived from single-molecule experiments

As shown above, the current results indicated the existence of an internal potential of the enzyme, a potential dependent only on the distance between Tyr639 and the palm domain, which facilitates translocation in real situations. To confirm the existence of this potential, the potential of mean force (PMF) was calculated directly during Tyr639 moved from PRE to POST positions in system 7, in which only the enzyme was present, using the umbrella sampling method. The calculated results are shown in Fig. 6(a). At the PRE position of Tyr639, which was located at d ≈ 0, PMF ≈−0.5 kcal/mol. At the POST position, which was located at d ≈ 0.4 nm, PMF ≈ −2.5 kcal/mol. The difference between the values of PMF at PRE and POST positions was about 2 kcal/mol . Fitting the available single-molecule experimental data to the general Brownian ratchet model of translocation resulted in a difference of about between the values of the free energy at pre- and post-translocation states with exclusion of the interaction between the enzyme and the DNA-RNA (supplemental information, text S2). The calculations of the difference in PMF were consistent with fitting the data to the free energy difference. This strongly supports our conclusion that the small free energy bias toward the post-translocation state is caused primarily from the difference between the internal potential energy of the enzyme itself in the pre-translocation state and in the post-translocation state. As shown in Fig. 6(a), in addition to the minimum value at the POST position, PMF also showed a minimum value at about 0.2 nm away from the POST or PRE position. This is also consistent with the conclusion that there exists an ITM for which Tyr639 is located between the PRE and POST positions.

Fig. 6. (color online) Free-energy profiles of translocation. (a) The potential of mean force (PMF) during the process of Tyr639 moving from PRE to POST positions in system 7 where only the enzyme is present. The x axis is the moving distance of Tyr639 from PRE to POST positions and the y axis is PMF. The PRE position corresponds to d = 0, with PMF = −0.5 kcal/mol, and the POST position corresponds to d = 0.4 nm, with PMF = −2.5 kcal/mol. (b) The potential of mean force (PMF) during the process of Tyr639 moving from PRE to POST positions in system 3 where the αβ methylene ATP and the template nucleotide i + 1 are removed in the insertion complex (PDB 1S76). The x axis is the moving distance of Tyr639 from PRE to POST positions and the y axis is PMF. The PRE position corresponds to d = 0, with PMF = −2 kcal/mol, and the POST position corresponds to d = 0.4 nm, with PMF = −6.5 kcal/mol.

PMF during the process of Tyr639 moving from PRE to POST positions in system 3, from which the αβ methylene ATP and the template nucleotide i + 1 had been removed from the insertion complex (PDB 1S76) was also calculated. The calculated results are shown in Fig. 6(b), where the most stable state of the system is located at POST position of Tyr639. At PRE position, which was located at d ≈ 0, PMF ≈ −2 kcal/mol. At POST position, which was located at d ≈ 0.4 nm, PMF ≈ −6.5 kcal/mol. The difference between the values of PMF at PRE and POST positions is about 4.5 kcal/mol . This PMF difference for system 3 is about 4kBT larger than that for system 7, which is due to the interaction between the side chain of residue Tyr639 and the i-th base pair. The ITM state can also be noted from Fig. 6(b).

3. Discussion

Comparison of the pre- and post-translocation structures indicates that the pivoting rotation of the O helix drives Tyr639 to move from the PRE to POST positions, inducing the translocation of the DNA-RNA relative to the enzyme by one nucleotide (or by a distance of 3.4 Å), which is called power-stroke mechanism.[8] Because the release of PPi triggers the pivoting rotation of the O helix, the power-stroke mechanism requires that the transition from the pre- and post-translocation states is located between the PPi-bound and PPi-release states. However, the results of fitting to the single-molecule experimental data obtained using the optical trapping method indicated that the translocation occurs after PPi release and before the next NTP binding.[18] The NTP binding acts as the “pawl”, biasing the rapid forward–backward translocation equilibrium between the pre- and post-translocation states, indicating that the Brownian ratchet mechanism was involved.[18] Fitting available single-molecule experimental data to the general Brownian ratchet model also produced a small free energy bias ΔE toward the post-translocation state rather than the pre-translocation state (supplemental information, text S2). In general, two models related to this small free energy bias ΔE have been put forward. The first proposes that the small free energy bias ΔE is caused by the pushing by the side chain of Tyr639.[19] The other proposes that the small free energy bias ΔE is due to the slightly larger binding affinity of the enzyme to the DNA in the post-translocation state than in the pre-translocation state.[31] However, the origin of the small free energy bias ΔE remains unknown.

Here, all-atom MD simulations were used to study the dynamic relationship between the movements of Tyr639 and the rotation of the O helix along with the rotation of the fingers domain. To make the simulation realizable, the (i + 1)-th base pair or the DNA-RNA substrate was removed. In this way, the movement of Tyr639 from PRE to POST positions becomes very fast, in the timescale of tens of nanoseconds. This simulation condition approximately corresponds to cases in which the large resistance to the movement, which results from the large binding affinity of the palm and thumb domains for the DNA-RNA substrate, is excluded. In real situations involving translocation, the palm and thumb domains moved relative to the DNA-RNA towards Tyr639, or Tyr639 and DNA-RNA moved towards the palm and thumb domains (Figs. 1(b) and 1(c)). All this occurs within a timescale of milliseconds. The movement requires overcoming the large binding affinity of the palm and thumb domains for the DNA-RNA. In this way, the mean translocation time in the real situation, , versus the mean movement time of Tyr639 from PRE to POST positions in these simulations, , is approximately , where ΔG is the large binding affinity of the palm and thumb domains for the DNA-RNA. For example, would give a about 105-fold longer than . It is here noted that, in the current simulations, the time required for Tyr639 to move from PRE to POST positions was significantly lower than that in the real situation. The effect of the rotation of the O helix on the movement of Tyr639 was not different than under real-world conditions. In this way, the simulations can be used to determine how the rotation of the O helix affects the movement of Tyr639 and how the movement of Tyr639 affects the rotation of the O helix and the fingers domain, elucidating the dynamic relationship between the translocation and the rotation of the O helix and fingers domain.

The current studies showed decisively that the PRE to POST movement of Tyr639 is not tightly coupled to the pivoting rotation of the O helix along with the outward rotation of the fingers domain. The two processes were found to be weakly dependent upon each other, suggesting that the power-stroke mechanism is not in play. Results also revealed the existence of a previously unknown ITM state for which Tyr639 was located about halfway between the PRE and POST positions. The existence of this ITM state has potential implications for the facilitation of NTP binding and selectivity. Moreover, results showed that an internal potential of the enzyme itself, dependent only on the relative distance between Tyr639 and the palm domain, produces the small free energy bias toward the post-translocation state relative to the pre-translocation state. The calculated value of the difference in internal potential energy was quantitatively supported by the available single-molecule experimental data. These findings suggest a modified Brownian ratchet translocation mechanism, as detailed below.

Following the completion of incorporation of NTP opposite the (i + 1)-th DNA base on the template, an elongation cycle begins with the active site in the palm domain of the enzyme positioned at the 3′-end of RNA primer. The O helix and the fingers domain are in the closed conformation and Tyr639 is in the PRE position (Fig. 1(b)). After the release of PPi, the O helix and fingers domain rotate rapidly to the semi-closed conformation. The RNA 3′-terminus occupies the active site, sterically blocking NTP from binding. Driven by thermal noise, the palm and thumb domains move forward (transitioning from the pre- to post-translocation states, with Tyr639 transitioning from the PRE to POST position) and vice versa (Brownian motion). The difference in internal potential energy ΔE of the enzyme facilitates the forward movement (Figs. 1(b) and 1(c)). Independently, the O helix and the fingers domain rotate to the open conformation. After the forward movement to the post-translocation state (Fig. 1(c)), the active site is left unoccupied, allowing the NTP to bind to the active site in the pre-insertion state. The occupation of the active site by NTP sterically prevents the palm and thumb domains from moving backward (the ratchet) (Fig. 1(d)). The transition from the pre-insertion to insertion states drives the N-terminal of the O helix to move towards the active site via the interaction between NTP and residue Arg627 in the N-terminal of the O helix (Fig. 1(a)). Simultaneously, the strong interaction between the side chain of Tyr639 and the newlygenerated base pair i + 1 drives Tyr639 outwards to the PRE position. This motion of Tyr639 (Fig. 1(d) to Fig. 1(a)) induces this increase in the internal potential energy of the enzyme by ΔE (note that this increase in energy is derived from the transition between the pre-insertion and insertion states). By comparison, the transition from the pre- to post-translocation states (Fig. 1(b) to Fig. 1(c)) induces the decrease in the internal potential energy of the enzyme by ΔE. After the incorporation of NTP into RNA, the next elongation cycle begins.

4. Materials and methods
4.1. Simulation methods

Three elongation complex (EC) crystal structures of the single-subunit T7 RNAP (PDB 1S76, 1S77, 1H38) were used for these MD simulations. The crystal structure 1S76 is the insertion complex and the crystal structure 1S77 is the pre-translocation complex, for which the fingers domain is in the closed conformation and Tyr639 is in the PRE position. The crystal structure 1H38 is the post-translocation complex, in which the fingers domain is in the open conformation and Tyr639 is in the POST position. All MD simulations were carried out by using GROMACS4.6[32] and the AMBER03 force field[33] was used.

The simulation methods and conditions used were the same as those in a previous work.[34] For convenience, they are briefly re-described as follows. All data regarding crystal structures used in these simulations were downloaded from the RCSB protein data bank. If there were any missing atoms, Swiss-Pdb Viewer was used to add them. A box was set to build a simulation system. To avoid the edge effect, the distance between the complex and the boundary of the box is at least 1 nm. To simulate the real system, solvent must be added, and necessary ions must be added with favorable concentration. Counter-ions were added to neutralize the system. All MD simulations were run at 300 K and 1 bar. The size of each step in the simulation was 2 fs and the neighbor list was updated every 5 steps. All chemical bonds were constrained using the LINCS algorithm.[35] The cutoff for van der Waals interaction and short range electrostatics interaction was 1 nm. The particle mesh ewald (PME) algorithm was used for calculations of long-range electrostatics.[36] Velocity-rescaling temperature coupling[37] and Berendsen pressure coupling[38] were used. The reference coordinates were scaled using the scaling matrix of the pressure coupling. Energy minimization was performed by using the steepest descent method for 50000 steps. The system was then equilibrated for 100 ps at 300 K and 1 bar pressure for 100 ps successively. In all MD simulations, residues 412–553 of the palm domain were restrained and all nucleic acids present in the elongation complex were left unrestrained. We have checked that without restraint of any residue, our simulations give similar results to those with residues 412–553 being restrained (Fig. S3).

The potential of mean force (PMF) was extracted from a series of umbrella sampling simulations.[39–41] A series of initial configurations was generated, each corresponding to a location wherein the residues of interest were harmonically restrained at increasing center-of-mass distance from a reference residue using an umbrella biasing potential. This restraint allowed each residue of interest to move within the configurational space in a defined region along a reaction coordinate between itself and its reference residue. To generate the initial configurations, all atoms of residue Tyr639 were pulled away from PRE to POST positions over 3 ns, using a spring constant of and a pulling rate of 0.2 nm/ns. Increasing the pulling rate 10-fold produced similar initial configurations, thereby ensuring that the pulling rate of 0.2 nm/ns would be sufficiently slow, as in previous works.[34] From the configurational trajectory, snapshots were taken to generate the starting umbrella sampling windows. The spacing of the simulation windows was less than 0.1 nm and the simulation period lasted 10 ns, so that the histograms of the configurations would sufficiently overlap with their neighbor windows. The simulation period of 10 ns is long enough to ensure the convergence of our results (see supplemental information, Fig. 4). Results were analyzed using the weighted histogram analysis method (WHAM).[42]

4.2. System establishment

Eight different systems were used for these simulations based on the high-resolution crystal structures. The systems were placed in cuboid boxes of simple point charge (SPC) water, in which there were 6 mM magnesium ions and neutralizing counter-ions. For system 1, the size of the box was 10 nm × 10.4 nm × 10.3 nm, with 30790 water molecules, 33 sodium ions, and 4 magnesium ions. For system 2, the size of the box was 11.6 nm × 10.9 nm × 10.2 nm, with 37991 water molecules, 25 sodium ions, and 5 magnesium ions. For systems 3, 4, 5, and 6, the size of the boxes was 10.2 nm × 10.4 nm × 10.3 nm, with 31603 water molecules, 31 sodium ions, and 4 magnesium ions. For systems 7 and 8, the size of the boxes was 10.2 nm × 10.4 nm × 10.3 nm, with 32173 water molecules, 12 chloride ions, and 4 magnesium ions.

Reference
[1] Steitz T A 2006 Embo J. 25 3458
[2] Landick R 1997 Cell 88 741
[3] Uptain S M Kane C M Chamberlin M J 1997 Annu. Rev. Biochem. 66 117
[4] von Hippel P H 1998 Science 281 660
[5] Steitz T A 1999 J. Biol. Chem. 274 17395
[6] Steitz T A 2009 Curr. Opin. Struc. Biol. 19 683
[7] Steitz T A Wang J Eom S H Jaeger J Restle T Jeruzalmi D 1996 Faseb J. 10 795
[8] Yin Y W Steitz T A 2004 Cell 116 393
[9] Kennedy W P Momand J R Yin Y W 2007 J. Mol. Biol. 370 256
[10] Tahirov T H Temiakov D Anikin M Patlan V McAllister W T Vassylyev D G Yokoyama S 2002 Nature 420 43
[11] Bar-Nahum G Epshtein V Ruckenstein A E Rafikov R Mustaev A Nudler E 2005 Cell 120 183
[12] Guajardo R Sousa R 1997 J. Mol. Biol. 265 8
[13] Guo Q Sousa R 2006 J. Mol. Biol. 358 241
[14] Vassylyev D G Artsimovitch I 2005 Cell 123 977
[15] Abbondanzieri E A Greenleaf W J Shaevitz J W Landick R Block S M 2005 Nature 438 460
[16] Bai L Santangelo T J Wang M D 2006 Annu. Rev. Bioph. Biom. 35 343
[17] Thomen P Lopez P J Heslot F 2005 Phys. Rev. Lett. 94 128102
[18] Thomen P Lopez P J Bockelmann U Guillerez J Dreyfus M Heslot F 2008 Biophys. J. 95 2423
[19] Yu J Oster G 2012 Biophys. J. 102 532
[20] Feig M Burton Z F 2010 Biophys. J. 99 2577
[21] Huang X H Wang D Weiss D R Bushnell D A Kornberg R D Levitt M 2010 Proc. Natl. Acad. Sci. USA 107 15745
[22] Kireeva M L Opron K Seibold S A Domecq C Cukier R I Coulombe B Kashlev M Burton Z F 2012 Bmc Biophys. 5 11
[23] Silva D A Weiss D R Avila F P Da L T Levitt M Wang D Huang X H 2014 Proc. Natl. Acad. Sci. USA 111 7665
[24] Woo H J Liu Y Sousa R 2008 Proteins 73 1021
[25] Golosov A A Warren J J Beese L S Karplus M 2010 Structure 18 83
[26] Radhakrishnan R Schlick T 2004 Proc. Natl. Acad. Sci. USA 101 5970
[27] Yang L J Beard W A Wilson S H Broyde S Schlick T 2002 J. Mol. Biol. 317 651
[28] Xie P 2008 Biosystems 93 199
[29] Pettersen E F Goddard T D Huang C C Couch G S Greenblatt D M Meng E C Ferrin T E 2004 J. Comput. Chem. 25 1605
[30] Meng E C Pettersen E F Couch G S Huang C C Ferrin T E 2006 Bmc Bioinformatics 7 339
[31] Xie P 2012 Proteins 80 2020
[32] Hess B Kutzner C van der Spoel D Lindahl E 2008 J. Chem. Theory Comput. 4 435
[33] Cornell W D Cieplak P Bayly C I Gould I R Merz K M Ferguson D M Spellmeyer D C Fox T Caldwell J W Kollman P A 1995 J. Am. Chem. Soc. 117 5179
[34] Duan Z W Xie P Li W Wang P Y 2012 Plos One 7 e36071
[35] Hess B Bekker H Berendsen H J C Fraaije J G E M 1997 J. Comput. Chem. 18 1463
[36] Darden T York D Pedersen L 1993 J. Chem. Phys. 98 10089
[37] Bussi G Donadio D Parrinello M 2007 J. Chem. Phys. 126 014101
[38] Berendsen H J C Postma J P M Vangunsteren W F Dinola A Haak J R 1984 J. Chem. Phys. 81 3684
[39] Patey G N Valleau J P 1973 Chem. Phys. Lett. 21 297
[40] Torrie G M Valleau J P 1974 Chem. Phys. Lett. 28 578
[41] Torrie G M Valleau J P 1977 J. Comput. Phys. 23 187
[42] Kumar S Bouzida D Swendsen R H Kollman P A Rosenberg J M 1992 J. Comput. Chem. 13 1011